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Abstract 

This paper is concerned with the accurate numerical approximation of the spectral proper- 
ties of the biharmonic operator on various domains in two dimensions. A number of analytic 
results concerning the eigenfunctions of this operator are summarized and their implications 
for numerical approximation are discussed. In particular, the asymptotic behaviour of the first 
eigenfunction is studied since it is known that this has an unbounded number of oscillations 
when approaching certain types of corner on domain boundaries. Recent computational results 
of Bjorstad and Tjostheim ||, using a highly accurate spectral Legendre-Galerkin method, have 
demonstrated that a number of these sign changes may be accurately computed on a square do- 
main provided sufficient care is taken with the numerical method. We demonstrate that similar 
accuracy is also achieved using an unstructured finite element solver which may be applied to 
problems on domains with arbitrary geometries. A number of results obtained from this mixed 
finite element approach are then presented for a variety of domains. These include a family of 
circular sector regions, for which the oscillatory behaviour is studied as a function of the inter- 
nal angle, and another family of (symmetric and non-convex) domains, for which the parity of 
the least eigenfunction is investigated. The paper not only verifies existing asymptotic theory, 
but also allows us to make a new conjecture concerning the eigenfunctions of the biharmonic 
operator. 



1 Introduction 

In recent years there has been growing interest in the spectral theory of higher order elliptic operators. 
However, in contrast to the well-understood theory for second order operators, the higher order theory 



can be quite different and is certainly less well developed. The recent survey article 0] contains 
an account of many of the key results. The purpose of this paper is to explore, using reliable 
and accurate numerical techniques, some of the properties of the eigenfunctions of the biharmonic 
operator on domains in 3ft 2 . In particular we investigate the existence of so-called nodal lines in 
the neighbourhood of certain corners of the domain and also how the parity of these eigenfunctions 
may change with domain geometry. As will become clear from the quantitative description of these 
problems below this is a demanding computational task since very high numerical accuracy is required 
in order to resolve the phenomena under investigation. 

The biharmonic eigenvalue problem typically comes in two different forms: the clamped plate 
eigenproblem 

A 2 u = Xu, m 
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and the buckling plate eigenproblem 

A 2 u = \Au , (2) 

each on a domain Q C 3? 2 subject to appropriate conditions on the boundary dQ. Throughout this 
paper we shall consider only the zero Dirichlet conditions 

Ou 

u = — — = Vx G dQ. (3) 

on 

In the remainder of this introductory section we describe details of the specific problems that are to 
be considered. Section |2| then contains an outline and justification of the numerical techniques, based 
upon the use of conforming C° finite elements, that are used to undertake these investigations. This is 
followed in Sections |3| and ^ respectively by our findings concerning oscillations of the eigenfunctions 
in the neighbourhood of corners of the domain and the dependence of the parity of the eigenfunctions 
on the domain geometry. The paper concludes with a brief discussion of our results along with a 
consideration of the issues associated with the validation of numerical simulations such as these, 
including mesh convergence and the use of rigorous enclosure techniques. 



1.1 Oscillatory properties of eigenfunctions 



For each of the two problems (|1|) and (fj) it is shown in |6[|,[l2"f that in the neighbourhood of a corner 
of dQ with sufficiently small internal angle 9 any eigenfunction changes sign an infinite number of 
times. Moreover the ratio of the distance from the corner, along its bisector, of consecutive zeros 
of the eigenfunction tends to a limit which is dependent upon 9. It transpires that this ratio is an 
increasing function of 9 and is such that the oscillations cease to be present when the internal angle 
exceeds some critical value 9 C . 

These results are special cases of more general theorems concerning higher order elliptic operators 



in greater than one dimension JT8[. Their derivation is based upon an asymptotic expansion of the 
eigenfunction, u, centred at, and in a neighbourhood of, the corner. The leading term of such an 
expansion as one approaches the corner along the bisector of the angle (which we always do) is of the 
form cr p , where r is the distance from the corner, and p is a solution of the transcendental equation 

sin((p + l W =[) 
sin 9 

Assuming that the coefficient c is non-zero the relevant solution is the one with smallest real part. 
From this equation, which has only real solutions for 9>9 c = 0.8128tt = 146°30', it may be concluded 
that when the internal angle is greater than 9 C the eigenfunction has no sign changes near the corner. 
Conversely as 9 — > both the real and imaginary parts of the exponent of the leading term in 
the asymptotic expansion grow unboundedly, implying that the eigenfunction oscillates with high 
frequency but is damped out quickly for small 9. Moreover, when we let s n be the distance along 
the bisector from the corner to the n th zero of u, where n increases with decreasing distance from 
the corner, then it follows immediately from the asymptotic expression that 

~ e 71 ^ as n — ► oo, (5) 



s n+l 



where p = a + i/3 is the solution of @ with smallest real part. Furthermore, when r n is the distance 
along the bisector to the n th extremum of u and t n is the magnitude of this extremum, it may also 
be shown that a 

$n i tn I $n \ /«\ 

~ and ~ (bj 
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as n —y oo. The first goal of this research is to verify this behaviour for a sequence of domains which 
take the form of sectors of a circle with increasing arc length. In view of the above discussion and 
the data in Table [3] below it is clear that such a verification will require very high accuracy from the 
numerical methods used. This table presents the solutions p = a + i/3 of (|) as a function of 9, as 
well as the ratios e n ^ and e an ^ which appear in (^) and (||). These solutions of (f|) are obtained by 
Newton's method with double precision arithmetic. 



9 


10° 


20° 


30° 


40° 


50° 


60° 


70° 


a = 5fte(p) 


25.141144 


13.079480 


9.062965 


7.057831 


5.857356 


5.059329 


4.491404 


j3 = 9m(p) 


12.864086 


6.384388 


4.202867 


3.095366 


2.416840 


1.952050 


1.608491 




1.27662 


1.63571 


2.11169 


2.75918 


3.66884 


4.99972 


7.05073 


e <W/3 


463.97239 


623.95258 


875.20459 


1291.07923 


2026.03589 


3437.12115 


6452.99420 




9 


80° 


90° 


100° 


110° 


120° 


130° 


140° 


a = Ke(p) 


4.067435 


3.739593 


3.479215 


3.268096 


3.094139 


2.949023 


2.826869 


(3 = $sm(p) 


1.339586 


1.119024 


0.930373 


0.762118 


0.604585 


0.446356 


0.261695 




10.43532 


16.56743 


29.27404 


61.69387 


180.5992 


1139.464 


163533.23 


e am/p 


13890.112 


36267.559 


126534.61 


709058.99 


9607097.6 


1033431725 


0.5472 • 10 15 



Table 1: The solutions of @ as a function of the angle 9. The ratios e 7 ^ and e a7r ^ are also given. 



1.2 Parity of eigenf unctions for non-convex domains 

Our other main goal is to investigate the dependence of the parity of the eigenf unctions of problems 
(EH) and (|2|) on the geometry of Q. Given a domain Q C 3? 2 which is symmetric about some line (which, 
without loss of generality, we will choose to be the y-axis), it is known that the least eigenfunction 
of the Laplace operator, subject to zero Dirichlet boundary conditions, is always positive and of 
multiplicity one | I3fl . Symmetry arguments now force it to be an even function of x. No equivalent 



result holds for the biharmonic operator. Indeed, it is the purpose of this work to demonstrate 
numerically that no such result holds for problems of the form (|I|) and (0) in arbitrary symmetric 
domains Q. This is achieved in Section || by considering a family of non-convex domains whose 
boundary is defined by the closed curves 

y = ±(c + x 2 -x 4 ) (7) 

for c > (see figure [I] below). Domain monotonicity arguments [[0| imply that the eigenvalues of the 



biharmonic operator increase as c — > and converge to the eigenvalues of the same operator on the 
region associated with c = 0. This region is disconnected and its eigenvalues have multiplicity two. 
Variational considerations now show that for small enough c > the low-lying eigenvalues occur in 
pairs, one associated with an even eigenfunction and the other associated with an odd eigenfunction. 
However in the biharmonic case one cannot say which of these is the smaller on purely variational 
grounds. Indeed by monitoring the smallest two eigenvalues of (|IJ) or (|2]) as functions of c, it may 
be observed that their paths cross on more than one occasion as c — > 0. Furthermore it is apparent 
that one of these eigenvalues corresponds to an even eigenfunction while the other corresponds to 
an odd one, and consequently we observe that the parity of the least eigenfunction is not fixed as a 
function of c. We again remark on the need for very high accuracy in our computational algorithms; 
this time in order to identify values of c for which the eigenvalue paths cross. 
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Figure 1: The dumb-bell shape domain for c = 1, c = 0.5, c = 0.1 



2 Numerical method 

Given the above discussion, it is clear that our discretization scheme must be able to accommodate 
quite general domain geometries. For this reason we choose a numerical scheme which is based upon 
the use of conforming C° Lagrange finite elements (see |21|, [p5|l for example) on some triangulation 



of the domain f2. The particular method that we use is based on that of Ciarlet and Raviart in |Tl 

Let H 1 ^) represent the usual Sobolev space of L 2 (Q) functions whose first partial derivatives 
are in L 2 (Q), and 



dn 



= . 

on ) 

Consider a family of regular and quasi-uniform triangulations, T h (where h is the diameter of the 
largest triangle), of the domain Q (see, for example, Ciarlet |10|), then we may define the space 

S h m = {p G C°(ti) : p\ T e P m {T), VT G T h }, 

where P m {T) is the space of polynomials of degree at most m over triangle T. From this we may define 
the following finite-dimensional subspaces of H 1 (Q): V h = and W h = D H^(Q). Then the 
Ciarlet-Raviart mixed finite element method approximates the solution u of (|J) by the component 
u h of the solution {v h ,u h } G V h x W h to the following problem. 

Problem 2.1 Find {v h ,u h } G V h x W h such that 

f v h Wl dtt+ ( Vu h ■ Vwt dVl = Vw! G v\ 

Jn Jn 

[ Vv h ■ Vw 2 dn = -X f u h w 2 dVl Vw 2 G w h . 
Jn Jn 



In this work we choose m = 3, i.e. piecewise cubic approximation, since such a discretization 
is known to satisfy the Brezzi stability condition, | 
property ( [[iT| ): 



p0| and has the following asymptotic error 



\u — u 



\m{n) + \\Au-v 



\L 2 (Q) 



< c\\u\\ H5{n) h 
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provided the solution u is sufficiently smooth. Hence the trial functions v and u h in Problem [2TT 

may be expressed as 



and 



u 



(9) 



where the </>.,■ 's are the usual piecewise cubic Lagrange basis functions on T h (see [pl],|25| for example), 
and rii and are the number of nodes in the interior of Q and on its Dirichlet boundary respectively. 
(We also use the convention that the interior nodes are numbered prior to those on the Dirichlet 
boundary.) This discretization leads to the following matrix problem 



M K l 
K O 



u 



-A 





Mu 



(10) 



where M G sftK+™d)x , M G $l n * xn * is a principal submatrix of M, K G sft"i x K+"d) anc j k* is 
the transpose of K. The specific entries of these matrices are given by 



Kij = J V0i ■ V0j dH, 



Mi. 



and the vectors v G Jj ni+n d and u G ^t ni denote the coefficients of v h and u h respectively in (|9]). 

It may be observed that u and A in ([TO]) are generalized eigenvectors and eigenvalues respectively 
for the corresponding Schur complement problem: 



KM- x K l u = XMu. 



(11] 



" M 


K* ' 











K 







u n+1 




-Mu n 



In this paper we are typically interested in approximating the smallest eigenvalues of this problem 
and so make use of straightforward inverse iteration (as described in |15 for example). Since the 
Schur complement matrix KM~ l K l is dense, whereas the matrix on the left-hand side of ( ]T0| ) is 
sparse, it is convenient to apply the inverse iteration to the latter system. This yields the following 
matrix equation at the (n + l) th inverse iteration: 



(12) 



When approximating the second smallest eigenpair, as in Section § for example, this iteration is 
slightly modified by subtracting out the component of the smallest eigenvector from u n+1 at each 
iteration. 

We solve the systems ( |T2"D using a direct sparse method based upon an initial block reduction 
followed by a sparse Cholesky factorization of symmetric positive definite sub-blocks as described in 
H. This permits re-use of the factorizations for the second and subsequent solves, which ensures 
that the cost of these is significantly less than that of the initial solution of (0). An additional 
advantage of a direct, rather than an iterative, method such as this is that it is possible to obtain 
highly accurate solutions: something that is required in both Sections [3] and i below. To demonstrate 
the accuracy of this approach (both the discretization and the direct solution), we now present some 
typical numerical computations for a biharmonic eigenvalue problem which has been considered in 
some detail elsewhere, fffll.H. 



Recall from Subsection |1.1| that in the neighbourhood of a corner with a sufficiently small internal 
angle, the eigenfunctions for both problems (|l|) and (§) are known to change sign at an infinite number 
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of points. In particular, when the internal angle 9 = tt/2 the solution of @ with smallest real part 
is 3.739593284 + 1.1190248i (see Table [I]). It then follows from (|) and (§) that the asymptotic ratio 
of distances along the bisector of consecutive zeros and extremal points of the eigenf unctions should 
both tend to 16.567428 as n —>■ oo. Also from (|(J) we see that the ratios between the magnitudes of 
consecutive extrema along the bisector tends to 36267.550 as n — > oo. In recent years many authors 
have tried to verify numerically these theoretical results using a variety of computational techniques 
(see, for example, P|,P|,|1T6|,|l24l|). Since the oscillatory features occur very close to the corners and 



are damped out very quickly, most of these attempts have, due to discretization errors and numerical 
inaccuracy, failed to find more than one sign change. A notable exception to this is the recent 
paper of Bj0rstad and Tj0stheim ([Q) in which the authors report five correct sign changes for the 
principal eigenfunction of problem (jl|) on the domain (0, 1) x (0, 1). This is achieved using a spectral 
Legendre-Galerkin method and by performing computations with quadruple precision arithmetic. 

In Tables ^| to [5] we present some comparative numerical results of our own, obtained using the 
mixed finite element method described above on three different meshes. In addition we contrast in 
the last column of the tables 0, [| and |5] the best results of Bj0rstad and Tj0stheim reported in ||. 
(This latter paper does not present results for the positions, s n , of the zeros and so no comparisons 
are possible for these figures). It may clearly be seen that, for each of these meshes, we are also able 
to identify at least five sign changes for the principal eigenfunction of ([D on the unit square. These 
calculations have been undertaken using double precision arithmetic on one eighth of the domain, 
using symmetries at x = 1/2, y = 1/2 and y = x, and with heavy local refinement of the meshes near 
the corner (h « 10~ 7 in the vicinity of the corner with a gradual transition to h »s 10~ 2 at the centre 
of the domain). The results are entirely consistent with the quadruple precision results presented in 
0] and with the predicted asymptotic ratios @ and (H). This leads us to have a reasonable degree 
of confidence in the underlying numerical procedure that we use for the investigations undertaken in 
the following two sections. 



3 Eigenfunctions on a circular sector 

The purpose of this section is to extend the numerical studies in which concern the behaviour 

of the principle eigenfunction of the biharmonic operator near the corners of a square domain, to the 
case of a general internal angle 9. This is achieved by considering a one parameter family of unit 
radius circular sector domains, Qg, and applying to (|1|) the finite element discretization defined in 



Problem |2.1| . Recall from Subsection that the oscillatory behaviour of the eigenfunction, u, near 
to a corner depends upon the internal angle 9 through the imaginary part of a solution of equation 
(|]). In particular, there is a critical angle, 9 C , above which no oscillations are present. As usual all 
of our calculations are restricted to the behaviour of the eigenfunction on the bisector of the angle 
at the corner. In this study we examine in detail the behaviour of the oscillations in the principal 
eigenfunction as 9 approaches 9 C = 0.81287T = 146°30', but we also present approximation of the 
eigenfunction for other values of the angle 9. 

It has already been demonstrated in Section ^| that our finite element discretization scheme is 
capable of producing the accuracy that is required to resolve up to five sign changes when 9 = n/2, 
on the unit square, which is consistent with the best numerical results of which we are aware 
Clearly a key to obtaining results with this level of resolution near to the corner is the use of an 
appropriate triangulation T h . For the circular sector regions that we consider here we use a similar 
mesh generation strategy to that adopted for the one-eighth of the unit square in the previous 
section. This involves the creation of an unstructured mesh based upon polar coordinates centred 
at the vertex. When r is greater than some (very small) chosen value, p\ say, the mesh is such that 
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FE method 


L-G method [4 


N 


5575 


22351 


90028 


5000 


Ai 


1294.934809196 


1294.934021432 


1294.933981245 


1294.933979592 


Sl 


0.042311747881 


0.042310932855 


0.042310963855 


— 


n 


0.032634299694 


0.032630435492 


0.032629530244 


0.032629472978 


h 


-0.169815505603- 10~ 4 


-0.169795463221 • 10~ 4 


-0.169791420686- 10~ 4 


-0.169791287981- 10~ 4 


S2 


0.002553961304 


0.002553862560 


0.002553860600 


— 


T2 


0.001969665622 


0.001969564641 


0.001969500077 


0.001969491936 


t 2 


0.468412256179- 10~ 9 


0.468173039472 • 10~ iJ 


0.468161662361 • 10~ 9 


0.468161006275- 10~ u 


S3 


0.000154151151 


0.000154149893 


0.000154149497 


— 


r3 


0.000118849965 


0.000118881375 


0.000118877347 


0.000118877352 


h 


-0.129136601165- 10~ 13 


-0.129088794006- 10~ 13 


-0.129085648295- 10" 13 


-0.129085369432- 10~ 13 


Si 


0.000009304900 


0.000009304402 


0.000009304373 






0.000007177212 


0.000007175502 


0.000007175314 


0.000007175365 


u 


0.356091995705- 10~ 1S 


0.355934162890- 10- 1S 


0.355925959506- 10- 18 


0.355925258182- 10" 18 


S5 


0.000000561678 


0.000000561628 


0.000000561569 






0.000000433175 


0.000000432837 


0.000000432811 


0.000000432767 


h 


-0.982438718508 • 10" 23 


-0.981777297429- 10^ 3 


-0.981303186480 -lO- 23 


-0.974860961611- 10" 23 



Table 2: Positions of the local zeros (s n ) and local extrema (r n ), and the values of the local extrema 
(t n ), for the first eigenfunction of the clamped plate problem QTJ) for n — 1 to 5, as calculated on 
three different meshes (where N is the number of unknowns in the system (|T0|) for each mesh). The 
last column contains the corresponding results from [4] where available. 



h ~ r/ii, and when h < p\ the mesh is approximately uniform. In this approach it is necessary to 
represent the circular arc by a piecewise affine approximation of side length In order to obtain 
an idea of the magnitude of the errors resulting from such an approximation we choose to perform all 
calculations on two different polygonal domains: one containing the sector and the other contained 
by it. It follows from the min-max principle that the m th eigenvalue of the differential operator on the 
sector is contained between the corresponding eigenvalue on each of these computational domains. 

In addition to considering the errors that result from taking an approximation to the domain Qg, 
it is also desirable to have an indication of the errors associated with approximating the continuous 
operator by the finite element discretization. In the numerical results that follow the calculations for 
each different value of 9 are presented for two different pairs of meshes of the type described above. 
Typically, the first pair of such meshes has only about 30% of the number of degrees of freedom as 
the second. Hence, when intervals for the local zeros/extrema positions obtained on the finer pair of 



JV 


5575 


22351 


90028 


Sl/s 2 


16.5671 


16.5674 


16.5675 


S2/S3 


16.5679 


16.5674 


16.5674 




16.5667 


16.5674 


16.5674 


s 4 /s 5 


16.5663 


16.5668 


16.5685 


limit 


16.5674 


16.5674 


16.5674 



Table 3: Ratios between consecutive local zeros as computed on three different meshes (where N is 
the number of unknowns in the system (|T0|) for each mesh). 
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FE method 


L-G method 4] 


N 


5575 


22351 


90028 


5000 


ri/r 2 


16.5684 


16.5673 


16.5674 


16.5675 


r2/n 


16.5727 


16.5675 


16.5675 


16.5674 


r 3 /r4 


16.5593 


16.5677 


16.5675 


16.5674 


r 4 /r 5 


16.5689 


16.5778 


16.5784 


16.5801 


limit 


16.5674 


16.5674 


16.5674 


16.5674 



Table 4: Ratios between consecutive local extremal positions as computed on three different meshes 



(where N is the number of unknowns in the system (10) and the system from [4] respectively). 





FE method 


L-G method g 


N 


5575 


22351 


90028 


5000 


\h/t 2 \ 


36253.4292 


36267.6722 


36267.6901 


36267.7125 


\h/t 3 \ 


36272.6177 


36267.5198 


36267.5223 


36267.5498 


\ts/U\ 


36264.9548 


36267.6044 


36267.5564 


36267.5496 


\U/t 5 \ 


36245.7209 


36254.0633 


36270.7433 


36510.3612 


limit 


36267.5596 


36267.5596 


36267.5596 


36267.5596 



Table 5: Ratios between consecutive local extremal values as computed on three different meshes 
(where N is the number of unknowns in the system ( |T0"D and the system from [4] respectively). 



meshes are contained within those obtained on the coarser pair, we can have a reasonable expectation 
of the reliability of these calculations. 

Tables ^| and |7| show the computed values of s n and r n (n = 1 to 4), the positions of the zeros 
and extrema respectively of the first eigenfunction of ([I]), for domains with 9 between 30° and 
140° inclusive. 

In each case the numerical results quoted take the form of pairs of values corresponding to those 
calculated on the inner (bottom value) and outer (top value) polygons respectively. It should be 
noted that, with one exception, each pair of values in Table |^ lies between the corresponding values 
in Table ||. While this proves nothing rigorous about bounds on the eigenvalue or the positions of 
the zeros and the extrema, it does provide a degree of confidence in the calculations. Furthermore, 
the one value (r 3 when 9 = 120°) for which there is any discrepancy between the two sets of results 
is clearly on the limit of our computational accuracy. Even in this case however, all four estimates 
of this value agree to two significant figures and have an absolute difference of less than 10~ 9 . The 
results in Tables || and |7| clearly show that, as the theory outlined above predicts, the frequency of the 
oscillations of the eigenfunction, u, decreases as 9 increases. Table [8] further illustrates this point by 
presenting the ratios of the positions of consecutive zeros and extrema for the fine mesh calculations. 
This table also shows the known asymptotic limit for these ratios in each case. It is notable that 
these asymptotic ratios are remarkably well approximated even by the first few ratios that we are 
able to calculate here. Furthermore, as 9 approaches an angle of 140°, it is clear that our numerical 
calculations are having great difficulty in resolving even the second change of sign. This is to be 
expected since the solution of equation (1) with smallest real part is equal to 2.826869 + 0.261695z for 
this value of 9: hence, by (|5]) and (|>D, the limiting ratio of the position and magnitude of consecutive 
extrema is equal to 163533.23 and 5.472 x 10 14 respectively. Given that r\ is already very small and 



S 



a 


30° 


40° 


50° 


60° 


70° 


80° 


hi 


0.065 


0.087 


0.110 


0.132 


0.154 


0.176 


Pi 


5 • 10- 4 


5 • 10~ 5 


1 • IO" 6 


2.5 • IO" 7 


1 • io-« 


1 • io- 8 


N 


34009 


32980 


36655 


33421 


34450 


30040 


Ai 


0^24.194986 
00 ^43. 169162 


i cq84. 340834 
i00 99. 470006 


oqOO.8711417 
° y 14. 1217006 


cr 7 16.2444537 
01 28.5030533 


4n 2i.2855619 
* u 33. 0287132 


q n 17.9563370 
ou 29.4734311 


Sl 


0.309733y28l 


231 n 79l3 ^ u 

u -^° iu 241412 


n 17n 76il453 
u - 1 ' u 6976522 


n 19 o7319839 
u - izo 6657365 


n nsfi9<J08096 
u.uou 8374821 


n n^s959seo5 

u.udoz 044100 


T\ 


n 9709386428 
U - z ' yZ 012640 


n onoi 650011 

u.zuzi 168919 


n i 4 r:2i9975i 


0.1024S 


0-07011S 


u.U4t>y 009493 


S2 


0.1459™ 


D flfiQK949688 

u.uooo 750755 


0-046™ 


0.0247 4 ^g^ 


0.0123™ 


n nn ^829245 

u.uuoo 776108 


T2 


1 31 fi54558Y 
u - i<5i0 369362 


0731 463995 

U.UMJ-289927 


0.0395™ 


0.0204™ 


0.0099^ 


n n n4 4024756 
U - UU4 3982853 


S3 


0.0691,™ 


0302 969386 


0.01267^36 


0.00494 894 * 9 


n noi 74798U7 

U.UUl '4 67069 


0.00053^ 




0.0623 2 2i67i 


0.02650§jgjjjj 


0.0107^ 


0.00409^$ 


n noi 41 14887 

U.UU141 04601 


0.000421 9 ^ 


Si 


0.03272^ 


0.0109™ 


0.00345^ 


0. 000989^| 


0.000247^ 


0.0000512^ 


r 4 


0.02951™ 


0.00960™ 


0.00293S 


0.000819^ 


0.000200^ 4 


0.000040^ 




e 


90° 


100° 


110° 


120° 


130° 


140° 


hi 


0.199 


0.222 


0.245 


0.268 


0.291 


0.315 


Pi 


1 • io- 8 


5 • 10~ 9 


2.5 • IO" 9 


5 • 10- 1U 


5 • 10- 1U 


5 • 10- 1U 


N 


29746 


24601 


23092 


22651 


20887 


19270 


Ai 


9 q76.3i55992 
z<5 87.7996658 


1Q 41. 3070303 
ly 52.8971428 


i fi 32.6947173 
J- 44 . 4979497 


1/L 05.7227746 
i4 17.8265343 


1 9 33.87885i3 
iz 46. 3582958 


1 100.6652772 
LL 13. 5880500 


si 


U.U^D4 024346 


n n2n4 y20930 


0.0096^^ 


003? V3V93U 
U.UU,3Z 6678 3 2 


0.00051SS5 


0.0000035^ 


n 


n n9SlU63243 
U.UZ5 072468 9 


0.0154S 4 


0071 498633 
U-UU/ 1369993 


0.00237^ 


0.00036$$? 


0.00000248^ 


S2 


0.00219^^ 


n n nn 7000095 

u.uuu 6989686 


0.000156jH$ 


0.000018 ^ 


0.00000045^ 




T2 


0.00169^ 


0.00052™ 


0.0001 15^ 


0.0000131^1 


0.00000032^ 




S3 


0.000132^ 4 


0.000023j$j$ 


0.00000253^ 


0.000000100^ 






r-i 


0.000102$;$ 


0.0000180^ 


0.00000187^ 


0.000000072^ 






S4 


0. 0000080 m 


0.00000081^ 


0.000000041J 








r 4 


0.0000061^ 


0.00000061^ 


0.000000030t 









Table 6: Local zeros (s n ) and local extrema (r n ), as calculated on the coarse meshes, of the first 
eigenvalue for the clamped plate eigenvalue problem on a circular sector of angle 9. 



\t\\ ~ 1.30729 ■ 10 -15 times the maximum magnitude of u, double precision calculations could not 
reasonably be expected to locate the second extremum with any significant figures of accuracy. For 
this reason, in Table || we are only able to present calculations of the location of the first zero of u 
when considering values of 9 greater than 140°. The parameters used for the generation of the meshes 
for these calculations are pi = 5 • 10~ 10 and hi ^0.16, which, in principal at least, allow us to observe 
a sign change at a distance of O(10~ n ) from the corner. In practice however we are unable to detect 
the first change in sign when 9 = 145° and 9 = 146° using such meshes. Nevertheless it is clear that 
these simple calculations do verify the theory which predicts the loss of oscillations at some critical 
angle 9 C , even though we are unable to accurately determine a value for 9 C using these meshes. 

Another interesting situation to investigate numerically would be the limiting behaviour of the 
principal eigenfunction as 9 — > 0. Theory predicts that, as 9 decreases, oscillations will become more 
frequent and their amplitude will decay more quickly as a function of r (although the asymptotic 
ratio of successive eigenvalues will decrease). This presents different, but equally demanding, com- 
putational challenges to those for 9 — > 9 C . In particular it is no longer sufficient to refine the mesh 
heavily only in the neighbourhood of the corner since changes in sign occur frequently throughout a 
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a 


30° 


40° 


50° 


60° 


70° 


80° 


hi 


0.033 


0.044 


0.055 


0.066 


0.077 


0.087 


Pi 


1 • io- 2 


2 • 10~ 3 


2.5 • 10" 4 


5 • 10~ 5 


1 • 10~ 5 


1 • 10~ e 


N 


82084 


82627 


87865 


87283 


86701 


90775 


Ai 


oca 28. 863922 
00 ^33. 606980 


i c;o88.056086 
i00 91. 837688 


o Qn 4.1158527 
oyu 7.4275436 


tr 7 19.2365745 
01 22.2999555 


4094-1419602 
^ uz 7.0760866 


qn90.7476335 
JUZ 3. 6247696 


Sl 


0.3097gftB 




0.1707293299 


n i 9 q7158333 
u - izo 6992733 


0.0868 8 ™ 


0.05 82326311 


n 


0.2792™ 


fl 9091466898 
u - zuzi 346638 


n 14 c;2ll6893 
U.140 lg81go6 


1 n24 451849 
U.1U/4 31471 3 


u - u ' ui 473169 


0.0459^ 


S2 


0.1459~ 


O.O83584935Q 


0.0465$$$ 


0.0247™ 


0.01232™ 


n055 8U1286 


T2 


0.13164^ 




0.03955$$ 


0.02048^ 


0.0099^ 


0.00440fS 


S3 


0.06911^ 


0.03029^ 


0.01267™ 


0.00494 82884 


0.001747^ 4 


0.000534 8 ,^ 4 


?~3 


0.06232™ 


0.02650S 


0.01078^ 


0.004097^ 


0.00141S 


0. 00042 1 8 ^ 


S 4 


0.03272$$$ 


0.01097™ 


0.003455^$ 


0.000989g^ 


0.0002478^ 


0.0000512^2 


U 


0.02951™ 


0.00960™ 


0.002938^ 


0.000819^ 


0.0002001^ 


0.0000404?$ 




9 


90° 


100° 


110° 


120° 


130° 


140° 


hi 


0.098 


0.110 


0.121 


0.132 


0.143 


0.154 


Pi 


2.5 • 10" 7 


5 • 10~ 8 


1 • io-« 


2.5 • 10" 9 


5-10- 10 


5 ■ 10- 10 


N 


88447 


87865 


87283 


85828 


85246 


79096 


Ai 


9 q79. 0885238 
^81. 9568311 


1 Q4 4.0948765 
iy4 6.9890139 


1 fi q5. 5227545 
L06 8. 4693656 


1 4 08. 6112832 
6320780 


19 q6. 8449411 
izc >9. 9585482 


1 1 n 3.7240407 
iiu 6.9471887 


Sl 


0.0364^^ 


0.0204™ 


0-0096S 


0.00327S 


0.000513^ 


0.00000353^ 


n 


0.02809S 


0.01546S 


0.00714^ 


0-00237S 


0.000366^ 


0.00000248^ 


S2 


0.00219gjj$g 


0.000699^^ 


0.000156 4 ^ 4 


0.0000181^ 


0.000000450^ 




T2 


0.00169S 


0.000528 4 j 2 | 


0.000115^ 


0.00001316^ 


0.00000032^ 




S3 


0.0001327^ 


0.000023^ 


0.00000253^ 


0.000000100^ 






^3 


0.0001023^ 


0.0000180^ 


0.00000187^ 


0.00000007^ 






S4 


0.00000801^ 


0.000000816^ 


0.0000000411 








r 4 


0.00000617^ 


0.000000616| 


0.0000000304 









Table 7: Local zeros (s n ) and local extrema {r n ), as calculated on the fine meshes, of the first 
eigenvalue for the clamped plate eigenvalue problem on a circular sector of angle 9. 



large part of the domain. Hence, for very small angles 9, a very dense finite element mesh is required 
over almost all of Q0 which, despite the narrowing of the domain, leads to very large algebraic systems 
(|10|). In fact, even for 9 = 10°, we observe a reduction in the apparent rate of convergence when 
considering the position of the zeros of u on sequences of finer meshes. This leads us to believe that 
a thorough numerical investigation of the behaviour of this eigenfunction as 9 — > will be an even 
more computationally demanding task. 

In the numerical results that we present above it has been possible, for the first time that we are 
aware, to verify a number of theoretical asymptotic results concerning the existence of nodal lines in 
the vicinity of a corner with arbitrary internal angle 9. Furthermore we have computed, with some 
confidence, estimates for the exact locations of zeros and extrema of the principal eigenfunction for 
a number of different values of 9. More comprehensive results, including those for the buckling plate 
problem, can also be found in |19[. In the following section we move on from providing numerical 
verification and extensions of known analytic results to demonstrate that numerics can also give some 
insight into problems for which there are no theoretical predictions. 
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e 


30° 


40° 


50° 


60° 


70° 


80° 


s\/s 2 


2.1219600 


2.7642773 


3.6711281 


5.0006316 


7.0510454 


10.435415 


S2/S3 


2.1121939 


2.7592690 


3.6688543 


4.9997170 


7.0507290 


10.435321 


S3/S4 


2.1116928 


2.7591817 


3.6688418 


4.9997176 


7.0507283 


10.435323 


limit 


2.1116888 


2.7591817 


3.6688419 


4.9997171 


7.0507288 


10.435321 


ri/r 2 


2.1210622 


2.7637792 


3.6708945 


5.0005086 


7.0510156 


10.435471 


r 2 /r 3 


2.1121602 


2.7592601 


3.6688530 


4.9997561 


7.0507030 


10.435314 


r 3 /r 4 


2.11168^ 


2.7591814 


3.6688418 


4.9997024 


7.0507288 


10.435311 


limit 


2.1116888 


2.7591817 


3.6688419 


4.9997171 


7.0507288 


10.435321 




e 


90° 


100° 


110° 


120° 


130° 


140° 


si/s 2 


16.567452 


29.274052 


61.693788 


180.59907 


1139.4^ 




S2/S3 


16.567430 


29.274040 


61.693978 


1 on c89y5 

iou.o 9586 






S3/S4 


16.567425 


29.274$| 


fi1 718979 
Di -673203 








limit 


16.56743 


29.27404 


61.69387 


180.5992 


1139.464 


163533.2 


n/r 2 


16.567440 


29.274017 


61.694352 


180.5981$ 


1139.4^ 




r 2 /r 3 


16.567535 


29.27394^ 


61.693^4 


180.62^ 








16.5674^ 


29.274*$ 


61 -751464 








limit 


16.56743 


29.27404 


61.69387 


180.5992 


1139.464 


163533.2 



Table 8: Ratios of consecutive local zeros and extrema of the eigenfunction (along with the theoretical 
asymptotic value) for different values of 9. 






141° 


142° 


143° 


144° 


145° 


146° 


147° 


Si 


0.00000111^ 


0.000000240$ 


0.0000000275^ 


0.000000000^ 









Table 9: Calculated positions of the local zeros of the eigenfunction for values of 9 above 140°. 



4 Parity property of eigenfunctions 

In this section we again consider problems (0) and (fj) but now on another family of domains, fl c say, 
which are defined as the regions bounded by the curves (0) for c > 0. We are particularly interested 
in the behaviour of the eigenfunctions which correspond to the smallest eigenvalues as c — > 0. In 
this limit the connected domain Q c tends to a disconnected domain which is symmetric about the 
y— axis. It follows therefore that each eigenvalue must have an even multiplicity and, in particular, the 
smallest eigenvalue will be repeated. When c > 0, but small, therefore, we expect the eigenvalues of 
the biharmonic operator to occur in distinct pairs (i.e. there will be a small gap between the {2k — l) th 
and the 2k th eigenvalues for positive integers k). The same argument holds when considering the 
Laplacian eigenvalue problem on Q c , subject to zero Dirichlet boundary conditions. Indeed, it is 
known that for this second order problem the eigenfunction corresponding to the lower eigenvalue in 
each pair is always an even function of x while the other eigenfunction in the pair is always odd. No 
such result is known in the case of the biharmonic operator however. In Table [10] below we present 
numerical evidence that such a result is false. This is achieved by demonstrating that branches of 
the two least eigenvalues (corresponding to eigenfunctions of different parity), viewed as functions of 
c, appear to cross on more than one occasion as c approaches zero. 
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c 


^&vg.7i / ^odd 






N 


i n 

J. . \J 


304005 


51 9393395 


170 8501 Q7 


1 29609 


n q 


349358 


68 981 4038 


1 q^ 448499 


141 1 99 


n 8 


41 0835 


q4 ^4377^9 


930 195849 


1 ^9649 


n 7 


D 4q4^q3 


1 3q 44D088 


981 q9874q 


1 08990 




608090 


991 940406 


365 091 690 


1 1 981 


5 


75501 Q 


385 377891 


51 491 666 


1 33634 


4 


914099 


794 064578 


799 107738 


112754 


n 35 


Q79477 


1 005 87754 


1 D34 34630 


1 90530 


n 395 


QQ050Q 


1 1 85 48005 


1 1 q6 83874 


1 98306 


3 


1 000849 


1 397 46767 


1 396 99991 


1 341 38 


n 975 


1 004906 


1 650 06659 

± \J O KJ . W U \J O £j 


1 649 01 01 7 


141914 


n 95 


1 nn4q43 


1 q55 47056 


1 945 851 90 


1 51 634 


n 995 

u — — ■ ' 


1 003936 


9330 q54qq 


9393 4371 5 


1 61 354 


n 9 


1 001465 


9800 48607 


9796 38901 


1 36514 


n 1 75 


1 000403 


3396 70646 


3395 33q41 


1 951 38 


n 1 5 


1 000031 


41 65 99544 


41 65 0q438 


1 36996 


0.125 


n qqqq«7 


5172 36322 


5172 42820 


152066 


0.1 


0.999998 


6518.06051 


6518.07172 


144002 


0.075 


1.000000 


8356.71703 


8356.71678 


136082 


0.05 


1.000000 


10936.0178 


10936.0178 


133634 


0.025 


1.000000 


14670.7338 


14670.7338 


126362 


0.01 


1.000000 


17752.3432 


17752.3432 


118082 



Table 10: Estimates of the two smallest eigenvalues, and their parity, for the clamped plate problem 
on Q c . 



The results presented have been calculated using the same piecewise cubic mixed finite element 
method described and used in the previous sections. In practice our computations were performed 
on a discretization of Q c in which the curved boundary is represented by an interpolating polygon. 
Since our motivation here is primarily to find an example in which the parity of the eigenf unctions 
changes with the domain geometry, this approximation is of no significant consequence. We again 
make use of unstructured meshes of triangles in the interior of the domain and use N to represent 
the total number of degrees of freedom in the finite element approximation. For each value of c 
calculations were performed on a sequence of meshes with different numbers of triangles. Only the 
results of computations on the finest meshes are presented in Table however an indication of 
the reliability of these figures is given by comparison with the corresponding calculations on coarser 
meshes. In particular, Table 11 gives the relative differences (as percentages) between consecutive 
ratios \ eV en/^odd computed on four different meshes for each value of the parameter c. In the table 
Ni (i = 1, ... ,4) denotes the total dimension of the linear system arising from the four different 
discretizations of the domain, and 6 (%) stands for the relative difference in the ratios \ e ven/^odd 
computed on consecutive meshes. It is reasonable to expect that our results are reliable if this 
relative "error" decreases (and becomes sufficiently small) as the mesh is refined. As we can see from 
Table [11] this is the case for all meshes and all values of the parameter c considered. 

Table ITU] clearly shows that, for our finite element discretization at least, there are indeed values 
of c for which the branches of the two smallest eigenvalues cross. Inspection of the corresponding 
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c 


1.0 


0.9 


0.8 


0.7 


0.6 


0.5 


0.4 


0.35 


0.325 


0.3 


0.275 


JVi 


52202 


54002 


59402 


41762 


47522 


53282 


50546 


54434 


57026 


60914 


63506 


N 2 


73442 


79922 


86402 


60482 


67394 


76034 


68042 


74090 


78626 


81650 


86186 


S(%) 


1.18c-2 


1.46e-2 


1.29e-2 


1.97e-2 


1.45e-2 


1.21c-2 


6.15e-3 


3.24c-3 


1.75e-3 


7.70e-4 


2.47e-4 


N 2 


73442 


79922 


86402 


60482 


67394 


76034 


68042 


74090 


78626 


81650 


86186 


N 3 


98282 


105842 


115922 


80642 


90722 


102818 


89858 


95042 


101954 


105410 


112322 


S(%) 


7.47e-3 


6.81e-3 


7.77e-3 


9.19e-3 


8.95e-3 


6.14e-3 


4.84e-3 


1.79e-3 


1.16e-3 


5.19e-4 


1.38e-4 


N 3 


98282 


105842 


115922 


80642 


90722 


102818 


89858 


95042 


101954 


105410 


112322 


N 4 


129602 


141122 


152642 


108290 


119810 


133634 


112754 


120530 


128306 


134138 


141914 


6(%) 


5.03e-3 


5.89c-3 


5.05e-3 


7.93e-3 


5.90e-3 


5.01e-3 


2.69e-3 


1.53c-3 


8.05e-4 


4.19e-4 


1.16e-4 




c 


0.25 


0.225 


0.2 


0.175 


0.15 


0.125 


0.1 


0.075 


0.05 


0.025 


0.01 




68690 


72758 


76466 


66818 


73730 


80642 


70562 


82658 


52562 


74882 


58466 


N 2 


92234 


96770 


95042 


84242 


93314 


102386 


92162 


108290 


75170 


91082 


76034 


S(%) 


3.02c-5 


8.81e-5 


4.88e-5 


2.42c-5 


4.80e-6 


9.99e-8 


9.99e-8 


0.000 


0.000 


0.000 


0.000 


N 2 


92234 


96770 


95042 


84242 


93314 


102386 


92162 


108290 


75170 


91082 


76034 


N 3 


119234 


127874 


115634 


103682 


112322 


126722 


116642 


122402 


101810 


108002 


95906 


S(%) 


2.01c-5 


6.43e-5 


3.64e-5 


1.79c-5 


3.10e-6 


9.99e-8 


9.99e-8 


0.000 


0.000 


0.000 


0.000 


N 3 


119234 


127874 


115634 


103682 


112322 


126722 


116642 


122402 


101810 


108002 


95906 


N 4 


151634 


161354 


136514 


125138 


136226 


152066 


144002 


136082 


133634 


126362 


118082 


S(%) 


1.41c-5 


4.42e-5 


2.78e-5 


1.50e-5 


2.60e-6 


0.000 


0.000 


0.000 


0.000 


0.000 


0.000 



Table 11: Percentage discrepancies S between the ratios Xeven/Kdd obtained on four different meshes 
for each domain Q c (discrepancies of less than 1.0e-8 have been rounded to zero). 



eigenf unctions demonstrates that one of these functions is always even and the other is always 
odd: we refer to the associated eigenvalues as A e?)en (c) and \ dd(c) respectively. Further numerical 



computations on a sequence of coarser grids, [|19|| , suggest that as the mesh size h is reduced, the 



locations of these crossing points remain stable and the difference between A ei)en (c) and X dd(c) does 
not tend to zero for values of c between these points. This suggests that it is reasonably safe for 
us to conclude that the parity of the eigenfunction corresponding to the least eigenvalue of (Q) does 
indeed change as the domain Q c evolves. Although we have only presented results for the clamped 
plate problem here, it is shown in |TS[ that similar behaviour may be observed for the buckling plate 



problem (|2]) on these domains. It is interesting to hypothesize on the behaviour of the eigenvalue 
branches as c gets very close to its limiting value of zero. The resolution of the computations that we 
have undertaken only appears to allow us to identify the first two or three values of c for which the 
least eigenvalue is repeated. Beyond this point discretization and rounding errors make it impossible 
to distinguish between them, even though they are almost certainly different for most values of c > 0. 
There may however be further values of c, smaller than those that we have been able to identify, for 
which \ e ven( c ) = ^odd(c): we conjecture that there are an infinite number of such values. 



5 Discussion 

In this paper our aim has been to demonstrate that the finite element method may be successfully 
applied to biharmonic eigenvalue problems of the form (Q) and ([J), and that it may be used both 
to verify existing asymptotic results and to allow new conjectures to be made about the behaviour 
of their eigenfunctions. An outline of the finite element implementation is given in Section ||], with 
complete details available in ||. Verification of the accuracy of this approach is provided by com- 
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parison with known analytical results applied on a unit square domain, along with the recent high 
precision numerical results (calculated using an entirely different discretization scheme) that appear 

in a. 

An advantage of the mixed finite element method used in this work is that it may easily be 
applied on unstructured grids of triangles, T h , and therefore on arbitrary polygonal domains in 3ft 2 . 
We have made use of this fact to investigate the spectral properties of the biharmonic operator on 
two further geometries: sectors of the unit circle and a family of non-convex domains. Furthermore 
it is possible to generate meshes for which the degrees of freedom in the corresponding piecewise 
cubic trial spaces are not uniformly distributed throughout Q. Where appropriate, this allows us to 
obtain a greater resolution of the eigenfunctions in those regions in which we wish to focus, without 
affecting either the dimension or the sparsity of the discrete linear system ([T0|). 

The results presented in Section ^| are in agreement with the asymptotic theory that appears in 



|6n,P2|,P7[|,P8| and also provide quantitative estimates of the positions of a number of sign changes 
and extremal points for various values of the internal angle 9. In Section f| we extend the numerical 
investigation further by considering a family of problems for which the eigenfunctions have signifi- 
cantly different qualitative properties. This has led to evidence that, unlike the Laplace operator, the 
parity of the least eigenfunction of the biharmonic operator can change as the domain is deformed. 

The unstructured finite element approach that we have used in this work may also be applied 
to further problems in this area. Clearly it is possible to investigate the spectral properties of the 
biharmonic operator on a greater variety of geometries in 2-d, and subject to a wider variety of 
boundary conditions. This would include, for example, further studies in the circular sector domains 
for very small angles 9. Perhaps more interestingly however (and certainly more computationally 
demanding) would be extensions to higher order operators and/or larger numbers of dimensions. 

We conclude this paper by noting that none of the computational results that are presented here 
provide rigorous information about the eigenfunctions considered (although the apparent mesh con- 
vergence observed is certainly a strong indicator). It would be interesting therefore to attempt to 
obtain provably correct bounds for certain quantities associated with the eigenvalues and eigenfunc- 



tions of the biharmonic operator on the regions considered. In ||24|| , for example, Weiners has made 
a start on this problem by computing enclosures for the first eigenfunction in the neighbourhood of 
a corner of the unit square. This involves computing estimates for the Sobolev embedding constants 
which in turn requires an upper bound for the defect of the eigenfunction. He is able to show that 
the defect is less than 0.000000893 which is sufficiently accurate to allow him to conclude that the 
first eigenfunction changes sign. Since the estimates of the amplitude of the second and higher sign 
changes that we have found are considerably smaller than for the first, new methods are likely to be 
required in order to compute the defect with sufficient precision. 
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